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ABSTRACT 


A proper understanding of the mechanism of contact between two or more 
nominally flat surfaces is crucial in the design process of many devices. This thesis, 
using analytical and computational methods, models the former through the use of fractal 
characteristics at the contact interface. A parametric analysis of the fractal surface was 
completed in order to properly understand fractal geometry and its effect on surface 
properties. The fractal surface was simplified so that Hertz theory could be used to 
model surface deformation and resulting contact stresses. The data gathered from the 
model was then input into an existing electromagnetic rail gun program to study the 
contact surface effect on exit velocity, temperature, electrical conductivity, and contact 
area ratio. Finally, a study of the fractal parameter effects on the electromagnetic rail gun 


was completed. 
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I. INTRODUCTION 


A. MOTIVATION 


In order to properly develop the technology used in the electromagnetic rail gun, 
one must take into consideration the importance of the contact mechanics involved. 
Mechanical, thermal, and electrical contact must all be understood at the interface 
between rail and armature. Figure | below displays a generic view of current rail gun 


construction. 





Figure 1. Rail Gun Structure [1]. 


The contact surface between the armature and the rails is an area that has been a 
problem for the development of the rail gun due to the simple fact that for proper 
operation, electrical contact must remain consistent throughout a firing. As the rail gun 
armature moves down the rails, the contact surface changes, and often results in changes 
in electrical contact. Issues arise when the surface profile changes enough so that an arc 
occurs. This arc causes a rapid increase in resistive heating well above the metals 


vaporization temperature which results in significant erosion [1]. 


This thesis investigates the mechanics of contact at the microscopic level so that a 
better understanding of the properties of nominally flat surfaces may be improved. Using 
the commercial computer program MATLAB, the mechanics of contact will be modeled 
and tested with the hopes of future improvements to the electromagnetic rail gun and 


similar defense-related issues. 


B. BACKGROUND AND LITERATURE SURVEY 


The first well-documented study of contact was that of Heinrich Rudolf Hertz in 
1881. The well known Hertz problem is one in which the contact stresses are solved for 
the simple cases of a sphere on a plane, a sphere on a sphere, and a cylinder on a 
cylinder; all perfect geometry [2]. Because Hertz theory is applicable for simple cases, 
many scientists in more recent history have developed more robust methods for solving 


contact problems. 


1. Greenwood and Williamson 


The goal of this research was to investigate the contact of nominally flat (but 
microscopically rough) surfaces. The assumption of nominally flat is defined as those 
surfaces in which the area of apparent contact is large so the individual contacts are 
dispersed and do not influence the neighboring contact locations. The purpose was to 
describe the theory of elastic contact and how it was closely related to real surfaces than 
the earlier theories. In order to compare the theoretical results, statistical results were 


tabulated and presented. 


A mathematical model was created to analyze the simple case of Hertz-contact (a 
spherical body in contact with a rigid flat plane) and determine the change in contact area 
with varying load. The mathematical model was then compared to experimental results 
to compare the change of pressure with different loads. Contact was then related to 
material hardness to compare the actual deformation with different loads. The 
development of a Taylor-Hobson model 3 Talysurf allowed the authors to analyze the 
surface topography of common materials. The research went on to describe the 


mechanics of material wear. 


Due to the fact that this research was completed in 1966, modern methods of 
calculation and computer-aided solvers were unavailable. The methodology began with 
the previous theory of Hertz contact and developed a surface model to analytically solve 


for deformation and area change with increased forces. Comparisons were made with 


experimental results. Also, statistical tabulation of the real surface profiles allowed the 
scientists to justify their assumptions of nominally flat surfaces as normal in most 


contact. 


The biggest conclusion of the research was that the separation between two 
surfaces depends on the nominal pressure (load divided by contact area), the number of 
micro-contacts and the total area of contact depend on the load only. These authors also 
developed the concept of “elastic hardness”, or basically that the area of contact can be 
predicted from the load the same way that it is for plastic contact using the conventional 
data for material hardness. They also discovered that most common surfaces have a 
Gaussian distribution of asperity heights [3]. The data gathered by these authors was 


cited as a reference for all of the other researchers mentioned in this thesis. 


2. J. F. Archard 


The purpose of this research was to investigate the hypothesis that elastic 
deformation of surface protuberances (asperities) was consistent with Amonton’s law, 
that the friction is proportional to the applied load. The basic idea investigated was that 
friction between sliding surfaces arises from local adhesion at the regions where 


protuberances come into real or molecular contact. 


Using Amonton’s law as a starting point, the author developed the theory of 
multiple small contacts in friction and backed up his results with experiments. The focus 
of the research was in the experimentation using a crossed-cylinders friction machine and 
the high polymer Perspex as the test material. Perspex was used due to its high capability 
to sustain large deformations under elastic conditions. The experimental results were 
used to develop coefficients of static friction based upon overall compressive load. 
Additional experiments were performed to analyze wear due to friction. The 
aforementioned experiments were redone using metal specimens and compared to the 
hypothesis that elastic deformation was the primary cause of change and not plastic 


deformation. 


Similar to Greenwood & Williamson, the author had limited computing power, 


and thus his work focused on experiments and tabulated results. Much of the research 
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focused on experimental setup, which required great precision in order to gather 
sufficient accurate results. The results of the experiments were compared to Amonton’s 


law of friction and results justified. 


The major conclusion from the research was that the real area of contact formed 
between two sliding surfaces was due to elastic deformation of the contacting 
protuberances. It was shown that the aforementioned assumption was consistent with 
Amonton’s law of friction. The conclusion was valid for both Perspex (with a large 
elastic range) and for carefully prepared metals under conditions where the damage was 


small [4]. 


3. Yan and Komvopoulos 


The goal of this research was to introduce a comprehensive contact mechanics 
analysis of elastic-plastic rough surfaces that can be characterized by three-dimensional 
fractal geometry and to present numerical results revealing the variation of the interfacial 


contact force and real contact area during quasistatic surface approach. 


The fractal surface equation was developed using the fractal geometry. Fractal 
geometry is common in nature, and can be observed in turbulence, precipitation, and 
specifically surface topography. Starting with the Weierstrass-Mandelbrot Equation 
(fractal equation), the engineers approximated the power spectrum, and generalized the 
equation using proper scaling properties. By including a random phase into the equation, 
they were able to develop a 3-dimensional function of x and y to create a surface with a 
fractal profile. Using this fractal profile, Hertz contact analysis was performed, with the 
major assumption that the contact areas could be modeled as individual spheres. Stress- 
Strain curves were used to model the elastic-plastic transformation. The numerical 
results were summarized to reveal the effect of material behavior on the resulting contact 


force. Graphs of the aforementioned numerical data were included. 


The methodology used in this research was a combination of 3 major areas. First, 
the overall approach was one that generated surface profiles using fractal mechanics. 
Using the fractal approach gave the authors the ability to randomly change the surface 


profile easily. The next major area was in the Hertz Contact Theory. This area, though 
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full of assumptions, was the focus of analyzing the actual contact between the two 
surfaces. Last, the results were compared statistically to ensure that the results were 


within the realm of actually observed experimental results. 


The research analyzed elastic-plastic rough surfaces characterized by 3- 
dimensional fractal geometry. The theoretical analysis of total contact force and real 
contact area in terms of surface separation distance, fractal parameters, and material 
properties helped to further develop insight into the effects of surface topography on 
contact. It was found that the actual contact area for elastic-plastic rough silica surfaces 
[the ones analyzed] was a remarkably small fraction of the apparent contact area, thus for 
most of the range of surface separation, the predominant deformation and resulting 


contact was due to elastic effects [5]. 


The research conducted by Yan and Komvopoulos was very much on the cutting 
edge of current contact research. The fractal approach was actually a major motivating 


factor for the research developed in this thesis. 


4. Sahoo and Ghosh 


The goal of this research was to use the finite element method to model the 
contact between 3-dimensional surfaces with fractal geometry. The focus was on 
developing non-dimensional contact area and displacement as functions of non- 
dimensional load for different surface profiles in the case of elastic contact. Also 


analyzed was the effect of strain-hardening in the case of elastic-plastic contact. 


Similar to the research of Yan and Komvopoulos, the Weierstrass-Mandelbrot 
function was used to model the surface topography. A finite element analysis was 
performed using the commercial software program ANSYS, which calculated the contact 
area after application of the load. A test of the finite element model was used with the 
simple sphere in contact with a rigid flat plane to ensure that the model could accurately 
determine the case of Hertz contact, which had already been completed experimentally. 
Both the analytical and FEM solutions were compared and solutions were generated for 


elastic contact and elastic-plastic contact. 


The methodology used by Sahoo and Ghosh was a combination of 4 major areas 
of research. First, the overall approach was one that generated surface profiles using 
fractal mechanics. Using the fractal approach gave the authors the ability to randomly 
change the surface profile easily. The next major area was in the Hertz Contact Theory. 
This area allowed the authors to develop an analytical solution. The third major area was 
in the use of the finite element method and commercial, off the shelf program [ANSYS] 
to develop a solution. The fourth, and possibly most important area of research, was the 
development of the solutions using non-dimensional methods which allowed their model 


to be independent size, scaling all numbers properly. 


The results of the FEM analysis agreed with the experimental results that were 
tabulated in previous literature. They discovered that in the elastic regime, contact area 
was linearly proportional to the contact load at small loads, but at high loads, the 
behavior became nonlinear. In the case of elastic-plastic contact, similar results were 
found. The system was linear for small loads, and nonlinear for larger loads with a direct 


effect on strain hardening [6]. 


This research is on the cutting edge of current contact analysis research. The 
fractal geometry matches nature closely, and can easily be scaled non-dimensionally. Of 
note by Sahoo and Ghosh, is that their solution is highly dependent on an accurate 
development of the true surface profile. Without knowing the specific details of the 
surface one is analyzing, the proper fractal surface can not be developed for accurate 


research. 


5. Ali and Sahoo 


The purpose of this research was to describe the adhesive contact between rough 
surfaces with small-scale surface asperities using an elastic-plastic model of contact 
deformation based on fictitious plastic asperity concept developed by other researchers, 
Abdo and Farhang [7]. The work began with a literature survey which discussed the 


work of Greenwood and Williamson and the downfalls of Hertz-contact assumptions. 


The topic of adhesion contact and the Tabor number (the ratio of elastic 
deformation to the range of action of adhesion forces) was introduced. The idea that the 
earlier analyses of adhesive contact treated everything on the macroscopic level as 
elastic-plastic, but at the microscopic level, the deformation of asperities was elastic then 
perfectly plastic. The concept of elastic-plastic adhesive forces was developed using the 
plastic asperity concept. Using a mathematical model, Ali and Sahoo were able to 
investigate three different asperity deformation conditions: purely elastic, elastic-plastic, 


and purely plastic. The results were solved and presented in graphical format. 


This modern approach to solving the elastic-plastic adhesive contact between two 
forces is at the forefront of using the computational method to solving difficult problems. 
Significant research was accomplished to investigate the shortcomings of prior research 
in order to develop a program that would solve for the details of different cases of 
deformation. Basically, where previous researchers assumed either plastic or elastic, Ali 


and Sahoo developed a program that could test the entire span of deformation. 


The results of the program were conclusive, and the prediction of load-separation 
behavior of contacting rough surfaces as functions of well-established elastic adhesion 
index and plasticity index. Basically, the results took previously tabulated material 
properties, which gave them the ability to analyze elastic and plastic behavior on both the 
macroscopic (material) and the microscopic (asperities on the surface). The 


representation of the data claimed to be a more-realistic model than all previous [8]. 


This research is the most current and up-to-date use of computing power to solve 
a problem that had previously been modeled. By taking known data for elastic or plastic 
contact and solving for the grey area in between, the authors created a program that is the 
most accurate representation of the real world. One suggestion that could be made for 
future work would be the inclusion of a graphical 3-D model in order to compare it to 


experimental deformation. 


6. Pratikakis 


The research of Nikolaos Pratikakis for his Master’s thesis was to develop a 


mathematical program to model the electromagnetic rail gun. His research used the finite 
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element method to solve for the potential field, temperature field, stress field, electrical 
conductivity, and thermal conductivity between the rails and the armature. His program 
used statistical data from experimentally gathered topographical data for typical rail and 
armature surfaces to model the contact interface. Using known analytical formulae for 
the multi-physics situation that develops within the rail gun, he was able to use the finite 
element method to determine the characteristics of the rail and armature during the entire 


cycle of a firing [9]. 


The Pratikakis program is used in this thesis to determine the properties of the rail 
gun when the armature and projectile exit the rails. His method was altered to account 
for fractal geometry at the interface between the rails and armature, but his finite element 
method solver was left as is. Using Pratikakis’ method, the thesis research conducted 


herein could be expanded for the application to the electromagnetic rail gun. 


Ts Polonsky and Keer 


The purpose of the research completed by Polonsky and Keer was to develop a 
faster way to solve rough contact problems. Arguing that the Hertz theory uses the 
stress-displacement response of individual microcontacts and does not take into account 
the effects of neighboring contact points, they sought to develop a method that would 


accurately describe generalized contact. 


The method that Polonsky and Keer adopted was the multi-level multi-summation 
and conjugate gradient techniques. They determined the real contact area using the 
conjugate gradient method. Then, using a two-dimensional multi-level multi-summation 
algorithm, they were able to keep the summation error under the discretization error. 
Their method views the contact surface as a three-dimensional grid in contact with 
another surface that has the same number of nodes on its respective grid. Assuming that 
the contact planes remain parallel to one another, their algorithm solves the three- 
dimensional general grid and then simplifies that grid into a two-dimensional grid that 
only represents the contact regions. Thus, their method is quick and accurate, taking into 


account the effects of regional contact points on one another [10]. 


This method takes into account regional contact effects on one another, but does 
so using purely elastic theory. It also models the contact as two dry surfaces in rough 
contact, which may or may not be the proper case in many lubricated mechanical devices. 
Further research and development would be helpful to develop a mathematical solver that 
takes into account lubrication, elastic-plastic deformation, and stress-strain relationships 


within the contact surface. 


8. Hu, Panagiotopoulos, Panagouli, Scherf, and Wriggers 


The objective of the research of Hu, Panagiotopoulos, Panagouli, Scherf, and 
Wriggers was to develop a program in which real contact areas could be found accurately 
using the finite element solution. Additionally, they sought to minimize computational 
costs while retaining accuracy. An adaptive finite element method was used to model 
fractal interfaces in a contact problem with the assumption of linear and finite elasticity. 
The focus on the research was in the problem description and the formulation of fractal 
interfaces. Contact kinematics was based upon different penetration characteristics when 
two different bodies came into contact with one another. This research concluded with 
the development and testing of the Finite Element Analysis Program (FEAP) using linear 


triangular elements [11]. 


The results and visual graphics of the developed program were similar to the 
commercial finite element program, ANSYS. Future development and research are 
recommended to improve the current program and allow for multiple types of contact to 


include elastic-plastic and friction. 


Cc, OBJECTIVE 


The objective of this research is to develop a mathematical model to determine 
the contact stresses between nominally flat surfaces. This situation is found in countless 
engineering applications, and is one of the prime factors in the development of the 
electromagnetic rail gun. By using simple Hertz theory applied to the complex geometry 
of a fractal surface, the intention is to create a simple model which will allow designers to 
choose different surface topographies and material properties during the development of 


the rail gun. The model developed will be an analytical model, strictly using 
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mathematics and computational analysis to determine contact stresses and area changes 
due to mechanical contact between two nominally flat surfaces. The commercial 
program MATLAB 7.2 will be used to develop all computer code. All assumptions will 


be clearly stated in appropriate sections. 
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YH. FRACTAL SURFACES 


A. FRACTAL REPRESENTATION 


The concept of “fractal dimension” was introduced by the mathematician Felix 
Hausdorff in 1918. This concept, distinct from the simple figures of Euclidean geometry 
was coined “fractal” by Benoit B. Mandelbrot of Poland. Mandelbrot was the first to 
point out the feasibility of modeling natural, physical objects with the concept. The 
property of self-similarity, which describes a system whose component parts represent 
the entirety of the system, make fractals an excellent choice for modeling surface 
topographies. Fractals remain invariant under changes of scale. The fractal concept has 
been adopted and used by mathematicians and engineers since 1975 [12]. The 
application of fractal geometry is applied hereafter to model physical topography with 


seemingly random features. 


B. CREATION OF BASELINE FRACTAL SURFACE 


The first step in modeling the fractal surface was to use the following fractal 
surface equation developed by Yan and Komvopoulos [5]: 


(D-2) V2 y n,, 
G Inv max (D-3)n 
= [| — oe > > 
(9) 2) [ M ” 


m=1 n=0 


ay" (x? _ yy” 


1 
7 (1) 


x 1COS@, , —COS 


To establish a baseline plot of the surface, MATLAB 7.2 was used with the 


following quantities substituted into the above equation: 





L = 1.0x10* (m) y = 15 ¢ = 5.9698 (rad) 
G = 136x110" (m) M = 10 (2) 
p= eg S000 








MATLAB code was developed to solve for the magnitude of z over a domain of 


one micrometer in the x and y directions, respectively. These heights, obtained by the 
equation for z(x, y) above, were then plotted using the embedded MATLAB code surf. 


This code creates a mesh plot with the spaces between the lines, called patches, filled in 


according to a color scheme based on the z data. 


In order to provide a cutaway view of the surface, the equation for z(x,y) was 


simply modified. The value of ” was held constant, arbitrarily at 50 nanometers, to 
display the profile in the middle of the sample. The values of < were then displayed 
using the MATLAB code plot. The results of these calculations are shown below as 
Figure 2. This two-dimensional profile more clearly shows the distribution of peaks and 
valleys in the domain. Of note is the randomness of the peak height distribution 


throughout the sample. 


MATLAB code was then developed to determine the location of the local peaks 
of the plotted surface. A peak was defined as a specific point whose value for Z was 
larger than both the preceding and following value of < in the x and y directions. The 
height of the peak, or < value, was then saved into a column vector for further 
calculations. Using this column vector, the peak radius was then solved for at each 
location deemed a peak. To solve for the peak radius, it was assumed that this radius was 
equal to the radius of curvature at the location of the peak. This was solved by setting the 


radius of curvature equal to the reciprocal of the second derivative as shown below: 


(Ax)’ 
+ x1) — 2% ay 


(3) 





r= 








(x41) 


This equation was solved in both the x and y directions, and then averaged to 


find the overall peak radius at that location. Although averaging in only the x and y 
directions may not be the most accurate assumption, it was used to simplify the 


calculations. 


12 


3-D Fractal Surface 
No. of Peaks = 1775 
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Figure 2. Baseline Fractal Surface. 


In order to gain a better understanding of the effects each parameter has on the 
surface profiles, the values for each of the variables in the z equation were varied. The 
following sections discuss the different parameters of the equation while displaying the 
results when each variable is increased and decreased from the baseline value. Included 
in the MATLAB plot were the mean peak height and mean peak radius over the domain 
shown. These values were calculated and displayed to assist in discerning between the 


different parameter changes. 
C. PARAMETRIC STUDY 


il Variation of the Number of Superposed Fractal Ridges, M 


The parameter M in the equation for z in terms of fractal notation represents the 
number of superposed fractal ridges. This parameter, if set equal to one, would construct 


a cylindrical, corrugated, two-dimensional surface. Thus, for values greater than one, the 
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parameter controls how many ridges are present in the surface. Figure 3 displays the 
fractal surface and profile when M is increased from the baseline value of 10 to 12. The 
change of M results in an increased total number of peaks over the domain, specifically 
from 1775 to 1835. The mean peak height increased, from approximately 7697 to 8601 
nanometers. Conversely, the mean peak radius is reduced from 0.59 to 0.18 nanometers. 
Thus, the peaks on average have increased height and decreased in radius, giving higher, 
narrower peaks. 
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Figure 3. Fractal Surface; Increased / . 


When the parameter M was decreased from the baseline value of 10 to 8, the 
equation yielded a surface and profile as shown in Figure 4. The decreased M resulted 
again in an increased total number of peaks from 1775 to 1822. The mean peak height 
decreased, from approximately 7697 to 6794. Similarly, the mean peak radius decreased 
from 0.59 to 0.19 nanometers. Thus the peaks on average have decreased height and 


decreased radius, giving lower, narrower peaks. 
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Although the three-dimensional surface profiles appear similar to the baseline, the 
two-dimensional profiles display the significant changes when the parameter M is 
changed. Increased values of M tend to reduce the difference between the peaks and 
valleys of the figure, causing the frequency band to become narrower. Surprisingly, a 
decreased value of M also caused some narrowing of the frequency band. The value of 
M was chosen as the baseline arbitrarily in order to generate a surface profile similar to 
one shown by Yan and Komvopoulos [5]. Because the equation was designed around the 
value of M equal to 10, its results yield a profile with the best distribution throughout a 


large frequency band. 


3-D Fractal Surface 
No. of Peaks = 1622 
Mean Peak Height =6794.647 (nm) 


2-D Fractal Profile at y = 50 (nm) 





Figure 4. Fractal Surface; Decreased M . 


Zs Variation of the Fractal Roughness, G 


The parameter G in the equation for z in terms of fractal notation represents a 
height scaling parameter. As it falls outside the summations and trigonometric functions 
in the equation, it is independent of frequency, and affects the results as a mere scale 


factor. It is termed the fractal roughness by Yan and Komvopoulos [5]. Increasing the 
15 


fractal roughness from the baseline 1.36x10"'' to 1.8x10'° meters simply shifts the 
surface and profile closer to the x- y plane. The overall number of peaks stays the same 
at 1775, but the mean peak height decreases from 7697 to 975 nanometers. Conversely, 
because the peak height is changing, it has an effect on the way that the peak radius is 
calculated. Thus, the mean peak radius increased from 0.59 to 4.69 nanometers, yielding 
shorter, wider peaks. The results of the increased value for fractal roughness are shown 


in Figure 5. 
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Figure 5. Fractal Surface; Increased G. 


Decreasing the fractal roughness from the baseline 1.36x10'' to 3.0x10°” 


meters simply shifts the surface and profile further away from the x-y plane. The 


overall number of peaks stays the same at 1775, but this time the mean peak height 


significantly increases from 7697 to 25791 nanometers. In addition, due to the 
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significant change in peak height, the mean peak radius decreased from 0.59 to 0.17 
nanometers, yielding shorter, wider peaks. The results of the decreased value for fractal 


roughness are shown in Figure 6. 


The parameter G, though it has an effect on the presentation of the surface 
profile, simply is just a scale factor. It keeps the same peak location, but varies the height 
and radius at that specific location. G is termed the fractal roughness because as it is 
increased, it creates lower, smoother peaks. The smoothness of the surface in this case, is 
the increased peak radius. The larger the peak radius, the greater surface area at the top 
of the peak; and thus an overall smoother appearance. 


3-D Fractal Surface 
No. of Peaks = 1775 
Mean Peak Height = 25791 .2493 (nm) 
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Figure 6. Fractal Surface; Decreased G. 


3. Variation of the Fractal Dimension, D 


The parameter D in the equation for z in terms of fractal notation represents the 
fractal dimension. Increasing the fractal dimension from the baseline 1.20 to 1.25 
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slightly increases the number of peaks from 1775 to 1780. The increased D decreases 
the mean peak height, however, from 7697 to 5640 nanometers. Conversely, it increases 
the mean peak radius from 0.59 to 0.87 nanometers. As D is increased, the topography 
is somewhat smoother, with a narrower frequency band and larger peak radii. Figure 7 


shows the new profile with D increased. 


When the fractal dimension is decreased from baseline 1.20 to 1.15, it has an 
opposite effect on the topography, yielding a much more jagged profile. The decreased 
D develops a mean peak height of 10510 nanometers, more than 2800 nanometers over 
the baseline. In this case the peak radius decreases from 0.59 to 0.39 nanometers. 


3-D Fractal Surface 
No. of Peaks = 1780 
Mean Peak Height = 5640.6154 (nm) 


height (nm) 


x (nim) 


2-D Fractal Profile at y = 50 (nm) 


height (men) 





Figure 7. Fractal Surface; Increased D. 


The variation of the fractal dimension, D, is a simple way to change the overall 
surface smoothness. By increasing D, the equation yields a smoother topography, with a 


narrower frequency band. Decreasing D has the opposite effect. The variation in fractal 
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dimension has little effect on the actual surface profile shape. The two dimensional 
profile remains nearly the same for the change in D. Specifically, the shape differences 
between the profile of the increased and decreased D values is negligible. Scale aside, 
the profiles are exactly the same. As shown in Figure 8, peak heights are in the same 
relative location and the overall number of peaks is roughly the same. 


3-D Fractal Surface 
No. of Peaks = 1769 
Mean Peak Height = 10510.3736 (nm) 


2-D Fractal Profile at y = 50 (nm) 





Figure 8. Fractal Surface; Decreased D 


4. Variation of the Density of the Profile Frequency, y 


The parameter y in the equation for z in terms of fractal notation represents the 
density of the profile frequency. Increasing the density of the profile frequency from the 
baseline 1.5 to 1.9 has no effect on the number of peaks which remains constant at 1775. 


The increased y very slightly decreases the mean peak height, however, from 7697 to 


7448 nanometers. It decreases the mean peak radius from 0.59 to 0.23 nanometers. As 


y is increased, the topography is somewhat smoother, with a narrower frequency band. 
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The peak height distribution is more consistent than in the baseline, meaning that many 
of the tall narrow peaks from the baseline have been replaced with shorter, narrower 
peaks when y is increased. Figure 9 shows the new profile with y increased. 
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Figure 9. Fractal Surface; Increased 7. 


When the density of the profile frequency is decreased, it has a greater effect on 
the profile. The number of peaks is decreased from 1775 to 1667. Decreasing the 
density of the frequency, however, causes the mean peak height to increase from 7697 to 
8520 nanometers and decrease the mean peak radius from 0.59 to 0.16 nanometers. 
Thus, the profile is shifted further from the axis, while at the same time the frequency 
band is decreased when compared to baseline. The surface, with increased profile 


frequency, is shown as Figure 10. 


Yan and Komvopoulos chose the value of y equal to 1.5 for their studies based 
upon prior research. They claim that this is the best value for the equation due to 
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considerations of surface flatness and frequency distribution density [5]. When compared 


to the increased and decreased y values, the baseline profile shows the most randomness 


and variation throughout the profile. 
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Figure 10. Fractal Surface; Decreased y. 
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Hl. CONTACT MODELING 


A. FORCE-DEFLECTION MODELING 


To properly understand the situation of contact, the problem was simplified so 
that simple Hertz theory could be applied. In order to apply Hertz theory, some 
assumptions were made in the solutions of the contact problem: 

e The contacting bodies are isotropic and elastic. 


e The contact areas are essentially flat and small relative to the radii of 
curvature of the undeformed bodies in the vicinity of the interface 


e The contacting bodies are perfectly smooth, and therefore only normal 
pressures need to be taken into account [2]. 


Taking into account these assumptions, well-known analytical Hertz equations 
were used to model each surface peak as a perfect sphere in contact with a rigid flat plate. 
This assumption required the development of a computer program to solve for each of the 


thousands of peaks in an extremely small cross-sectional area. 


iL: Hertz Equations 


The simplest form of the Hertz problem is one that describes two spherical 
surfaces in contact. The contact area in this case is a perfect circle and is described as 


having a radius, a, as follows [2]: 





P(E,+E,)nr | 
+ LY, 
«088 eee | (4) 
EE, G; + r,) 
The maximum contact pressure, o-,, is described as: 
P 
C,=la— > (5) 
qa 


Also described is the relative displacement, 6, of the centers of the two spheres 


due to local deformation and contact at the interface as shown below: 
4 1/3 
6 =0.77| P* : + : ye (6) 
E, E, ae) 
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When the aforementioned equations are again simplified to represent a sphere in 
contact with a rigid flat plate, the following substitutions are made: 
BE, =E,=E 
her (7) 
fr, =o 


Thus, the equations are simplified as shown below: 


1/3 2 \1/3 » \i3 
a=038{ 27") : a, = 06n| . sais) (8) 





r 2E*r 
For this analysis, however, the requirement to impose a deflection and measure 
the force required to cause that deflection was desired. Therefore, the deflection equation 


above was solved for the force, P , giving: 


3/2 


P =0.74{E\(dr'”) (9) 


Using the above equation for the radius, a, of the contact circle, the individual 


contact area for each sphere against a flat plate was solved as follows: 


3]? 
A=2a’ anx{oss(22 


Using this equation and the MATLAB program, a force-deflection analysis was 


completed. 


2. MATLAB Program 


The simplicity of the Hertz theory allows the development of an analytical 
solution to be relatively straight forward. Using a complex surface generated by the 
fractal equation described in the previous chapter, the program analyzes a specific area 
and identifies the peaks. The peaks are searched for in both the x and y directions 
individually and are tabulated in a column vector based upon their respective coordinate 
value. Each time a peak is found, a number is added to the counter, which displays the 


total number of peaks in a given sample area. 
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At each of these locations deemed a peak, the height is recorded as the value in 
the z direction of the standard Cartesian coordinate system. A similar column vector is 


created to tabulate the local radius of curvature in the x and y directions respectively. 


These radii of curvature of each peak are averaged as a single radius of a sphere for 
simplicity. In order to use Hertz theory in the simplified single sphere on a flat plate, the 
peaks are analyzed individually and treated as local spheres with radii equal to the 


respective radius of curvature. 


The force required to induce (or resist) a given deflection was the desired goal of 
the program. To achieve this, the relative deflection needed to be solved at each of the 
peaks. But, since the overall surface had multiple peaks with multiple heights and 
different radii, a loop needed to be created to separate rigid body motion from contact 


deflection. 


The separation between the fractal surface and the rigid flat plate was defined as 
the difference in height between each individual peak and the highest peak on the surface. 
To clarify, if a rigid flat plate was suspended above the surface with no contact, no force 
would be required. This would be the case as the plate was lowered until it made contact 
with the highest peak. From this point on, an opposition force would be generated as the 
plate was lowered and more and more peaks came in contact. Thus, the highest peak was 


chosen as the reference for the separation data. 


To create a force-deflection plot over a span of increasing deflection, a loop was 
generated to step-increase the deflection from zero to an arbitrary maximum. Simply 
stated, the program simulated the action of pushing a rigid flat plate onto the peaks of the 
fractal surface, all the while remaining parallel to the x-y plane. To account for the 
different peak heights, the separation was used to determine when a specific peak would 
come into contact with the rigid flat plane. Some of the shorter peaks would never come 
into contact with the rigid plane, and thus were discounted. The program included a 
counter to determine the number of effective peaks; the number of peaks that actually 


were affected by the contact with the rigid plane. 
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Based upon the separation, contact, and changing relative deflection, calculations 
were completed at each individual effective peak. The opposition force at each peak was 
tabulated in a column vector and then summed to determine the overall effective force of 
opposition through the use of a loop. The summation data was tabulated in a column 
vector based upon each step decrease of 6, so that the force-deflection data could be 
plotted Additionally, using the equation for the contact area of each peak, the program 
determined the overall effective contact area for each step of 6. Using the data from the 
force and area calculations, a local contact pressure was solved at each peak and 
tabulated in a column vector to calculate the overall opposition pressure as the deflection 


was changed. 


The program described above was then slightly modified to include an extra loop 
that allowed the user to change some of the parameters of the fractal equation and plot the 
changes. This visual display of the different parametric changes helped to discover some 
direct relationships between surface profile shape and the changes in effective contact 
area and opposition force. In all the following figures, the zero displacement starts from 


the highest peak of the fractal surface. 


3. Results, Deflection-Generated Opposition Pressure 


Using the aforementioned MATLAB program, plots of the deflection-generated 
opposition pressure were created to determine the effects of each individual parameter. 
The figures on the following pages show the results when different parameters are 


increased and decreased from the baseline. 


a. Variation of the Number of Superposed Fractal Ridges, M 


As previously discussed in the parametric study, the number of superposed 
fractal ridges, M , has a significance on the shape of the fractal surface. An increase in 
M caused an average increase in peak height, and a decrease in peak radius. The 
expected effect on the deflection-generated opposition pressure matched the results as 
shown in Figure 11. Increasing the number of superposed fractal ridges caused a 


significant increase in opposition pressure. This increase is believed to have come from 
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an increased number of peaks, thus more peaks in contact. More peaks in contact 
resulted in a greater pressure distribution. Of note, however, is that along with the effect 
of the wider frequency band of the increased M , the larger deflections have a much 


greater effect on the contact pressure. 
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Figure 11. Pressure-Deflection Curves, Changing Parameter ™ . 


The taller peaks make contact first, but as the deflection is increased in 
steps, the pressure increases as more and more peaks become in contact. The slope of the 


pressure-deflection curve is much steeper for the case of the increased M . 


Conversely, when M is decreased, the effect on the fractal surface is one 
which decreases the frequency band, peak height, and peak radius. The decreased width 
of the frequency band has the most significant effect. Because many of the peaks are 
shorter and closer in height to the surrounding peaks, they all make contact at smaller 


deflections. Instead of changing rapidly as in the case of the increased M , the slope of 
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the pressure deflection curve remains much flatter, with very small changes in slope 


along the length of the curve. The effect of the decrease in M is shown in Figure 11. 


b. Variation of the Fractal Roughness, G 


The fractal roughness, G, simply represents a height-scaling parameter, 
and thus the effect of an increase or decrease is directly correlated to the pressure- 


deflection curve shown in Figure 12. 
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Figure 12. Pressure-Deflection Curves, Changing Parameter G. 


The Increasing the fractal roughness shifts the pressure-deflection curve 
up. The value for the peak radius is calculated based upon a change in peak height in a 
specific direction, and since an increase in G causes the peaks to become shorter, the 
radius is directly affected. When G is increased, peak radius increases. Wider, shorter 


peaks will have a greater contact area and thus greater pressure for the same force. 
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Conversely, a decrease in G has the opposite effect. Peak heights are 
significantly increased, but radius is significantly decreased. The taller, narrower peaks 
make contact in the same manner as the shorter wider peaks, but the decreased radius 
causes a decrease in contact area. Thus, the contact pressure is decreased, yielding a 


much flatter pressure-deflection curve as shown in Figure 12. 


Cc. Variation of the Fractal Dimension, D 


The fractal dimension, as previously discussed, changes the smoothness of 
the fractal surface. An increase in D develops a smoother topography with narrower 
frequency band. In addition, the peak radii are increased, which causes more peak 
surface to be in contact. The small changes of D are shown in Figure 13. Although the 
curves are slightly different, they all follow the same trend. As previously discussed, the 
surface profile shape remains virtually the same for different values of D, the largest 
change being the peak height. An increase in D caused a decrease in peak height, and a 
decrease caused a significant increase in height. Thus, the shorter peaks made contact 
over a larger area, causing a steeper pressure-deflection curve, whereas the taller peaks 


had the opposite effect. 
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Figure 13. Pressure-Deflection Curves, Changing Parameter D. 


d. Variation of the Density of the Profile Frequency, y 





The density of the profile frequency had different effects on the surface 


topography as was previously discussed. Thus, the profile was inconsistent for different 


values of vy. Predicting the shape of the profile based on the y value was inconclusive. 


Thus, the originators of the fractal surface equation, Yan and Komvopoulos, chose the 


value equal to the baseline of 1.5 based upon their prior research. The pressure- 


deflection curves are shown below in Figure 14, but direct correlation to the effects due 


to the parameter y was not possible. 
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Figure 14. Pressure-Deflection Curves, Changing Parameter y. 


4. Results, Deflection-Generated Contact Area Ratio 
a. Variation of the Number of Superposed Fractal Ridges, M 


The effect of the various values of M was fairly simple to determine 
based upon the results of the parametric study previously conducted. Increasing the 
number of superposed fractal ridges caused increased peak height and decreased radius. 
Although the initial peak radius had an effect on the contact area, the formulation of the 
equation to solve for the continued increase in contact radius after a peak made contact 
caused the actual area to increase significantly. The effect was that the taller peaks, made 


contact early, and their small peak radii increased at a rate much greater than the baseline. 


The results are shown below in Figure 15. 
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Conversely, the decrease in M caused the peaks to become shorter while 
at the same time having a decreased peak radius. Thus, the peaks came into contact much 
sooner and the peak radii would all increase at the same time as the other peaks in 
contact. Therefore, the percentage of area in contact did not significantly change over the 
range of deflections shown. The results matched what was expected based upon the 


fractal surface shape and profile. 


Deflection-generated Contact Area Ratio 
6=0 +50 nm 


0.25 


0.2 


0.15 


0.1 


% Astrea in Contact 


0.05 


° 
a 
weet" 
we 
eee” 


weneere® 





revartrs 


0 50 ~ 100 150 200 20 300 350 400 460 = 400 
Deflection (nm) 











Figure 15. Area Ratio-Deflection Curves, Changing Parameter ™ . 


b. Variation of the Fractal Roughness, G 


Similar to the effects upon the pressure-deflection curves, the variation of 
G had the exact same effect on the area-deflection curve shown below as Figure 16. The 
decreased G flattened the curve, whereas the increased G caused a steep increase in the 


curve. Again, the reason for this change in shape is due to the fact that the fractal 
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roughness simply is a height scaling parameter. The radius was calculated based upon 
the change in height (or radius of curvature) at the peak, and thus the effect is directly 


correlated. 
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Figure 16. Area Ratio-Deflection Curves, Changing Parameter G. 


Cc. Variation of the Fractal Dimension, D 


An increase in D developed a smoother topography with narrower 
frequency band. The smoother the topography is, the more area is in contact. When D 
is decreased, the topography becomes more jagged, and thus fewer peaks make contact. 
The decreased number of peaks in contact causes a decrease in contact area ratio. The 


results of the variation of fractal dimension are shown below as Figure 17. 
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Figure 17. Area Ratio-Deflection Curves, Changing Parameter D. 


d. Variation of the Density of the Profile Frequency, y 


Varying the density of the profile frequency, as previously stated, was 
inconclusive as far as estimation of surface shape and profile was considered. Thus, the 
effect on contact area ratio was also varied based upon the different profile frequencies. 
Again, the value of 7 equal to 1.2 was chosen by the originators of the equation based 
upon previous research in surface flatness and frequency distribution. The results of the 


area-deflection plot are shown below as Figure 18. 
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Figure 18. Area Ratio-Deflection Curves, Changing Parameter 7. 


5, Curve Fitting 


In order to properly analyze the results, it was desired to find suitable equations 
for the deflection-generated opposition pressure and the area-deflection curves, 
respectively. Using MATLAB and its embedded code for polynomial curve-fitting, the 
pressure and area curves were solved using linear, quadratic, and cubic equations. The 
results of each polynomial fitting were assembled and plotted along with the raw data to 
determine which equation best fit the data. The goal was to determine an equation that 
was both simple and accurate, and thus research beyond the level of the cubic polynomial 
was not considered. The results are shown below as Figure 19. Although the quadratic 
equation does a reasonably good job of fitting the data, it was determined that the cubic 
remains closer to the data throughout the range of deflection. These cubic equations were 


then used in the Railgun program discussed in the next chapter. 
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Figure 19. Polynomial Curve-Fitting. 
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V. RAIL-GUN MODELING AND ANALYSIS 


A. INTRODUCTION 


The research completed for this thesis focused on the contact problem 
encountered during the development of the electromagnetic rail gun. A study of the 
specific interface between the armature and the rails was necessary to properly 
understand the mechanics in this area. The thesis research of Nikolaos Pratikakis in 2006 
focused on the development of a mathematical model of the electromagnetic rail gun. 
Using the equations developed in the previous chapter, the Pratikakis model was slightly 


modified to reflect the research of this thesis. 


B. PRATIKAKIS MODEL 


The mathematical model created by Nikolaos Pratikakis was developed using 
MATLAB code and the finite element method. It was made up of two main parts which 


are briefly summarized below. 


i: Part I, Setup of Inertial Coordinate System and Geometric 
Configuration 


Part I of the program was the development of the specific geometric dimensions 
and mesh of the model based upon lengths and widths of the rails and projectile as 
provided by the manufacturer. The geometry was simplified into a two-dimensional, 
uniform rectangular mesh for both the rails and the armature. Of note, however, is that 


the geometry could be changed very simply to account for different shaped armatures [9]. 


2. Part II, Main Program and Description of Functions 


A schematic of the rail gun program is shown as Figure 20. This program 
completes many calculations during each loop that properly account for the multi-physics 
nature of rail gun functioning. Specifically, electrical current density and current 
distribution are used to determine the total heat generated due to electrical resistance. 
These values are also used to determine the Lorentz Force on the projectile. Using partial 
differential equations, the program also completes a stress analysis to estimate the x and 


37 


y displacements. Thermal stresses are also calculated by assuming that the rails are 
fixed and rigid. The acceleration of the armature is also calculated using the second law 
of Newton. The program continues to run in an iterative loop until the displacement of 


the projectile is larger than the length of the rails, at which it records an exit velocity [9]. 
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Figure 20. Schematic of Rail Gun Program [8]. 
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C. RESULTS OF MODIFIED PRATIKAKIS EQUATION 


In order to properly understand the effects of each of the fractal parameters in a 
real-world application, the Pratikakis rail gun program was used. Using this program, a 
baseline set of results was generated using the nominal values previously discussed in 
this thesis. Then, the following parameters were altered: number of superposed fractal 


ridges, fractal roughness, fractal dimension, and profile frequency. 


i: Variation of the Number of Superposed Fractal Ridges, 4 


Using 10 as the baseline value of M , the parameter was decreased to the values 8 
and 9, and increased to the values 11 and 12. Each time the parameter was changed, the 
Pratikakis rail gun program generated the values of the area ratio, electrical conductivity, 
maximum temperature, and velocity of the projectile at the time at which it exits the rails. 
These results were recorded for each of the four different variations of the parameter, M . 
In order to gain a better sense of the change, the results were all normalized against the 
baseline value for comparison. The results of this parametric change are shown as Figure 
21 with one exception. When the program was run with a value of M equal to 9, it 
produced results which were a significant outlier. When the normalized data was plotted, 
the M =9 data resulted in a plot that was almost impossible to discern the other values 


of M , and thus this set of results was left out of the figure. 


When the number of superposed fractal ridges was decreased or increased, it 
caused the values for all four of the measured output quantities to increase. As 
previously discussed, the decrease in M results in shorter, narrower peaks, whereas the 
increase in M results in taller, narrower peaks. The effect of the narrowing of peak 
radius may result in a greater contact area, which allows for improved electrical 
conductivity. The greater the electrical conductivity, the more efficient Lorentz force, 
which results in a greater exit velocity and higher temperature. Of note in Figure 21 is 
that all four plots display the same trend on different scales. This is believed to be a trend 


strongly influenced by the effective area ratio. 
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Figure 21. Variation of the Number of Superposed Fractal Ridges, M . 


Zz Variation of the Fractal Roughness, G 


Using 13.6 picometers (1.36 x10''m) as the baseline value of G , the parameter 
was decreased to the values 0.30x10"'' and 0.70x10"", and increased to the values 
1.50x10""' and 1.80x10° meters. Each time the parameter was changed, the Pratikakis 
rail gun program generated the values of the area ratio, electrical conductivity, maximum 
temperature, and velocity of the projectile at the time at which it exits the rails. These 
results were recorded for each of the four different variations of the parameter, G. In 
order to gain a better sense of the change, the results were all normalized against the 


baseline value for comparison. The results of this parametric change are shown in Figure 
22. 
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Figure 22. Variation of the Fractal Roughness, G . 


When the fractal roughness was decreased the area ratio, electrical conductivity, 
and temperature all decreased, but the exit velocity increased slightly at G=0.70x10"" 
and decreased at G=0.30x10"''. As previously discussed, when the parameter G_ is 
decreased, the surface peaks become slightly shorter and wider, but the number of peaks 
remains the same. The wider peak radii decrease the deflection that occurs at the contact 
surface, and thus decreases the change in contact area. This, in turn, decreases the 
conductivity, temperature, and exit velocity. An exception occurs at the value of 
G =0.70x10"'' meters, where the exit velocity increases when compared to baseline 


when all of the other parameters have decreased. 
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When the fractal roughness was increased, the area ratio, electrical conductivity, 
temperature, and exit velocity all increased. When G is increased, the result is taller, 
narrower peaks. Again, the effect of decreased peak radius is that the area ratio may 
change easily as an force is applied. The increased area ratio thus dominates the rest of 
the results, and yields an increase in electrical conductivity, temperature, and exit 


velocity. Of note, is when G is significantly increased as is the case when 
G =1.80x10° meters, resulting area ratio increase is immense, and has a direct effect 


on the other quantities calculated. 


3: Variation of the Fractal Dimension, D 


Using 1.20 as the baseline value of D, the parameter was decreased to the values 
1.15 and 1.175, and increased to the values 1.225 and 1.25. Each time the parameter was 
changed, the Pratikakis rail gun program generated the values of the area ratio, electrical 
conductivity, maximum temperature, and velocity of the projectile at the time at which it 
exits the rails. These results were recorded for each of the four different variations of the 
parameter, D. In order to gain a better sense of the change, the results were all 
normalized against the baseline value for comparison. The results of this parametric 


change are shown as Figure 23. 


As was previously discussed, the fractal dimension is the simplest way to change 
the overall surface smoothness. Increasing D yields a smoother topography with a 
narrower frequency band. The relative shape of the profile remains unchanged, such that 
peaks remain in the same location, but their heights and radii may change. Thus, one can 
expect different results when the frequency band is narrowed. Because the fractal 
dimension predominantly effects peak height, and in turn, frequency band, the results of 


the Pratikakis program are difficult to determine based on this parameter alone. 


In the case of D=1.15 and D=1.25 all results decrease compared to the 
baseline, whereas in the other two cases tested, the results increase. The conclusions 
gathered from this result is that the peak height and frequency band may have 
counteractive effects on the contact area. In some instances, the peak height may cause 
fewer peaks to be in contact, but be easily deformable, and in other cases such as a slight 
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increase, their may be more peaks in contact, but be more difficult to deform. Thus, one 
must examine the size of the frequency band and the actual profile before trying to 


determine the direct effects of the variation of the fractal dimension, D. 
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4. Variation of the Density of the Profile Frequency, y 


Using 1.5 as the baseline value of y, the parameter was decreased to the values 
1.3 and 1.4, and increased to the values 1.7 and 1.9. Each time the parameter was 
changed, the Pratikakis rail gun program generated the values of the area ratio, electrical 
conductivity, maximum temperature, and velocity of the projectile at the time at which it 


exits the rails. These results were recorded for each of the four different variations of the 


43 


parameter, vy. In order to gain a better sense of the change, the results were all 


normalized against the baseline value for comparison. The results of this parametric 


change are shown in Figure 24. 
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Figure 24. Variation of the Density of the Profile Frequency, 7. 


The density of the profile frequency is a fractal parameter that has inconsistent 
results when the parameter is increased and decreased. In some instances, a decrease in 
y results in a smoother, narrower frequency band, whereas in others, it yields a more 
jagged, wider frequency band. Increasing the profile frequency density has similar 


effects. Nonetheless, for the specific values of v tested, the baseline profile had the most 
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randomness and variation throughout the profile and was chosen by Yan and 


Komvopoulos to be 1.5 for the best random fractal profile. 


Of note in Figure 24 is that the shape of the curve for area ratio and electrical 
conductivity follow similar trends. Although a different shape than the area ratio and 
electrical conductivity curve, the temperature and exit velocity curves follow a similar 
trend as well. The unpredictable nature of the parameter y has a direct effect on the 


similarly unpredictable results of the Pratikakis equation. 
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VI. SUMMARY AND CONCLUSIONS 


In this thesis, a computer program was developed to model nominally flat contact 
surfaces throughout the mechanism of contact. The use of analytical methods and Hertz 
contact theory allowed the program to solve for results analytically. Ultimately, a 
seemingly random profile is reduced to a finite number of Hertz contact problems that 
can be solved independently. By understanding each parameter of the fractal surface 
generator, one can tailor the surface to any known surface profile. Because many 
surfaces found in nature reveal surface profiles similar to fractal geometry, the 


application of this program has limitless potential. 


A parametric study was conducted in order to further understand the fractal 
surface equation. A baseline was chosen based upon previous work by Yan and 
Komvopoulos, and the following parameters were varied: the number of superposed 
fractal ridges, fractal roughness, fractal dimension, and density of the profile frequency. 
The effects of changing each individual parameter were compared to the baseline and 


plotted. 


The number of superposed fractal ridges had an effect on peak height, peak 
radius, and frequency band changes. Variation of the fractal roughness was determined 
to be a scale factor. The fractal dimension parameter had the greatest effect on peak 


height changes. The density of the profile frequency had inconclusive results. 


Use of a preexisting electromagnetic rail gun program allowed application of the 
nominally flat contact surface results. The Pratikakis program clearly displayed the 
effects of surface profile upon performance of the electromagnetic rail gun. Each 
parameter was studied to determine its effects upon not only the surface, but the 
application to a real world problem. The parameters that were changed when analyzing 
the effect on profile shape were changed in the same manner when running the rail gun 
program, specifically increasing and decreasing the parameter from the baseline. The 
effect of each parameter had a direct effect on surface shape, which in turn changed the 


results of the rail gun program. 
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VII. RECOMMENDATIONS 


This research provides a method to model contact surfaces using a fractal 
equation. It simplifies the complicated problem of seemingly random contact down into 
a simple Hertz analysis. Although the problem is calculated on the three-dimensional 
level, it may be inaccurate in the fact that each contact peak is modeled independently. 
The effect that one peak has upon another may cause significant changes in the fractal 
surface and the way in which it reacts. Future studies using the finite element method 
would provide more accurate results. Also, a thorough validation of the model against 


any experimental test data would be beneficial to refine the program. 
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